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ABSTRACT 

Here we present CAFein, a new computational tool for investigating radiative dissipation of dynamic 
tides in close binaries and of non-adiabatic, non-radial stellar oscillations in isolated stars in the linear 
regime. For the latter, CAFein computes the non-adiabatic eigenfrequencies and eigenfunctions of 
detailed stellar models. The code is based on the so-called Riccati method, a numerical algorithm 
that has been successfully applied to a variety of stellar pulsators, and which doesn't suffer of the major 
drawbacks of commonly used shooting and relaxation schemes. Here we present an extension of the 
Riccati method to investigate dynamic tides in close binaries. We demonstrate CAFein's capabilities 
as a stellar pulsation code both in the adiabatic and non-adiabatic regime, by reproducing previously 
published eigenfrequencies of a polytrope, and by successfully identifying the unstable modes of a 
stellar model in the /3 Cephei/SPB region of the Hertzsprung- Russell diagram. Finally, we verify 
CAFein's behavior in the dynamic tides regime by investigating the effects of dynamic tides on the 
eigenfunctions and orbital and spin evolution of massive Main Sequence stars in eccentric binaries. 
The plethora of asteroseismic data provided by the NASA's Kepler satellite, some of which include 
the direct detection of tidally excited stellar oscillations, make CAFein quite timely. Furthermore, 
the increasing number of observed short-period detached double white dwarfs (WD) and the observed 
orbital decay in the tightest of such binaries open up a new possibility of investigating WD interiors 
through the effects of tides on their orbital evolution. 

Subject headings: (stars:) binaries: general - stars: interiors - stars: oscillations - dynamic tides 



1. INTRODUCTION 

The current state and evolution of binary systems is 
affected by a wide range of physical processes, the un- 
derstanding of which is important in interpreting obser- 
vations. In this work we focus on dynamic tides, the 
tidal regime in which the free oscillation modes (eigen- 
modes) of one of the binary components can be excited by 
the companion's periodic tidal potential, with the driv- 
ing frequency being comparable to the stellar eigenfre- 
quencies. For this purpose, we have developed CAFein 
(Code for non- Adiabatic, non-radial Forced stEllar os- 
cillatloNs), a novel computational tool to investigate in 
detail the impact of dissipation of dynamic tides in close 
binaries. The efficiency of non-adiabatic, dynamic tides 
in exchanging angular momentum between the binary or- 
bit and the component spins in a binary depends on the 
tidal energy dissipation mechanism and on its strength. 
In addition, this new dynamic-tides tool can be used for 
the study of non-forced stellar pulsations. 

Thanks to NASA Kepler^s unprecedented photomet- 
ric accuracy, the effects of dynamic t ides have become 
readily visible in electr omagnetic data (jWelsh et al.ll201ll : 
[Thompson et al.ll2012( ). A phenomenal example among 
Kepler^s Objects of Interest (KOI) is KOI 54, a highly 
eccentric binary hosting two A stars; their light curves 
clearly reveal that the free oscillation modes of one or 
both stars are tidally excited (Welsh et al..,2011f ). The 
theoretical modeling of such features not only allows to 
further constrain the stellar and binary properties, but 
also potentially provides a probe to the stellar interiors, 
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otherwise hidden to direct electromagnet ic observations 
(IBurkart et al.ll2012bl: iRiller fc Lail[2Q12a . 

Pioneering investigations targeting the observed circu- 
larization periods of massive (O, B, F) binaries identi- 
fied radiative damping as the main mecha nism f or the 
dissipation of dynamic tides. In particular, IZahnI (|1975[ ) 
was the first to invoke radiative damping of tidally ex- 
citated 5-modes in these massive Main Sequence (MS) 
binaries. A systematic comparison between the circu- 
larization periods predicted by Zahn's theory and the 
observations showed that many systems circularized well 
above the theoretically predicted period, showing evi- 
dence for a more efficient tidal dissipation mechanism 
(iGiu ricin et all [l98l iNorth fc Zah3 iSool iMazeh et al.1 
2006). A highly pro mising explanation to th is discrep- 
ancy was provided by iWitte fc Savoniid ()2001f ) . The au- 
thors followed simultaneously the exchange of angular 
momentum between the stellar spins and the orbit due 
to tides and the evolution of the star's eigenfrequencies 
due to natural stellar evolution. This detailed analysis 
demonstrated that a binary can be locked into a resonant 
state for a prolonged period of time (the so-called reso- 
nance locking). Such long- lasting resonances can dra- 
matically speed up dissipation and hence tidal evolution 
of a binary's orbit, yielding better agreement with the 
observed circularization periods. 

Beyond the extensive w ork on tides in non-degenerate 
systems (see iZahnI '2008' for a review), recent studies 
have focused on the effects of dynamic tides on the or- 
bital evolution of detached binaries hosting white dwarfs 
(WD) and, in particular, double WDs (DWD). These bi- 
naries are widely recognized as important gravitational 
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wave (GW) sources: they are the most numerous and 
guaranteed sources for the next-generation of space- 
based detectors sensitive to low-frequency GW s (e.g. 
LISA, iDanzmann fc the LISA studv teamlll996l and ref- 
erenc es therein, and cLISA/NGO, Amaro-Seoan e et al. 
20121 see also ^cleman s_et_aL 20.01a.b, 2004 ^Liu et al. 
20101: iRuiter et al1l2010l for theoretical predictions), and 
are currently observed electromagnetically. In the past 
few years, th e Extr emely Low Mass WD su rvey ( ELM, 
IBrown et all IMol l20Tl IKihc et al.l [20Tfll Mm suc- 
cessfully quintupled the number of known detached 
DWDs expected to merge within a Hubble time, bring- 
ing the number of systems to 24 (see Table 4 by 
IKihc et al.|[20ll for a summary of the currently known 
systems) and discovering DWDs with periods down to 
- 12min (SPSS J 65133.33-F284423.3 hereafter J0651, 
iBrown et aH 120111 : iHermes et all 120121 ). The exciting 
J0651 system harbors a tidally deformed He WD eclips- 
ing a C/0 WD in a detache d binary. Since its o rbital de- 
cay was recently measured (IHermes et al.ll20l"2l ). J0651 is 
one of the cleanest astrophysical laboratories to test our 
understanding of tidal dissipation in these sources, and 
WD interiors. 

Recent work has focused on the adiabatic tidal 
excitation of free oscillation modes in C/0 WDs. 
iRathore et al.l ()2005D used polytropic models to repre- 
sent C/0 WDs in eccentric binaries and focused on the 
/-mode. They found that adiabatic tides can drive the 
modes to high amplitudes potentially beco ming nonlin- 
ear. Similar conclusions were reached by iFuller fc Lail 
()2011|) . who targeted adiabatic, dynamic tides in cir- 
cular binaries, using detailed C/0 WD models and fo- 
cusing on the excitation of modes in the close pas- 
sage through a resonance. These authors also found 
that mode excitations can cause significant deviations 
in the orbital evolution of DWDs from the pure point- 
mass assumption and are very important in the spin 
synchronization process. Their analysis, being limited 
to the adiabatic regime, does not include any dissipa- 
tion, which can limit the non-linear growth of the reso- 
nant modes. As a follow-up on the violation of l inearity 
found in the adiabatic treatment, IFuller fc Lail (|2012bl ) 
considered the tidal excitation of gravity waves in C/0 
DWDs treating dissipation via the so-called "outgoing 
wave boundary condition" (BC). Such a BC implicitly 
assumes the waves are damped at the WD's surface via 
radiative damping or non-linear effects; as a result the 
formation of standin g waves is p r event ed. They ob- 
tain results si milar tolFuUer fc Lail ()20 111) . Using a sim- 
ilar approach, IFuller fc Lail (1201231) extended their in- 
vestigation to He WDs, focusing mainly on the effect 
of tidal heating and its observational signatures. The 
authors found that tidal heat is likely deposited in the 
outer layers of the WD and that it can dominate the 
WD's luminosity for the shortest orbital period binaries 
(5. 1 5 min) . Recent inve stigatio ns by iValsecchi et al.l 
(poU ) and iBurkart et all (|2012aD are more focused on 
the effect of dyn amic tides on th e orbit al and spin evo- 
lution of DWDs. IValsecchi efaTI (|2012[ ) applied CAFein 
to a He WD model representative of the He component 
in J0651 to investigate the effect of linear, dissipative 
(non-adiabatic) dynam ic tides on its orbital evolution. 
IValsecchi et al.l (|2012D calculated the full tidal response 
of the WD as multiple modes are excited simultaneously 



for a w ide range of driving frequencies. IBurkart et al.l 
(|2012al) studied the effects of linear and non-linear dy- 
namic tides in DWDs with circular orbits hosting both 
He and C/0 WDs. In this study the dynamical tide is 
approximated as a superposition of standing waves and 
the WD response is treated as a simple harmonic oscil- 
lator with driving and dissipation. In the linear regime, 
the damping processes considered are thermal diffusion 
and turbulent convection, while in the non-linear regime 
the damping time for traveling waves is set by the g- 
mode group travel times. The radiative damping rate is 
also approximated in the quasi-adiabatic limit (from the 
star's adiabatic eigenf unctions, instead of the full non- 
adiaba tic eigenfunctions). We refer to IValsecchi et al.l 
(|2012f ) for a discussion and compa rison with the results 
presented bv IBurkart et al.1 ()2012al ). 

We have developed CAFein, a novel code to compute 
both non-adiabatic and non-radial stellar oscillations in 
isolated stars and forced stellar oscillations in close bi- 
naries. Following our understanding of tidal dissipation 
in non-degenerate stars, which is now able to explain 
the observed circularization periods in open cluster bi- 
naries (Zahn 1975; Wittc fc Savoniic 2001), we consider 
radiative damping to be the main mechanism to dissipate 
dynamic tides in stars with radiative envelopes. 

In this paper we describe in detail the mathematical 
and numerical implementation of CAFein and we present 
comparisons to past results found in the literature. In 
§ 121 we introduce CAFein as a stellar oscillation code. 
In § 12.1! we briefly summarize the equations governing 
non-adiabatic and non-radial stellar oscillations and in 
ii l2.2l we give a detailed description of the Riccati method 
implemented to solve such equations. In § 12.31 we test 
CAFein's results and accuracy by calculating eigenfre- 
qucncies and eigenfunctions of different stellar models. 
In § 12.3.1! we calculate the adiabatic eigenfrequencies of 
a polytrope and check how the results change if the rel- 
evant parameters entering the Riccati method are var- 
ied. In § 12.3.2! we compare these eigenferquencies with 
previously published results. In 12.3.31 we verify the re- 
liability of CAFein in the non-adiabatic regime by iden- 
tifying the unstable modes of a stellar model in the (3 
Cephei/SPB regiorQof the Hertzsprung-Russell (HR) di- 
agram. In § |3| we move onto dynamic tides and their 
implementation in CAFein. In § !3.1! we introduce the 
tide-generating potential. In § !3.2! we summarize the 
equations governing tidally excited stellar pulsations and 
the secular evolution of the orbital elements and stellar 
spin. In § !3.3! we describe the modifications applied to 
the Riccati method to solve the stellar pulsation equa- 
tions when the tide-generating potential is included, and 
test such extension i n §13.4! In § 13.4.1! we reproduce the 
results presented bv lPolfliet fc Smeverd ([1990) and show 
that a dynamical tide can be approximated as the sum of 
the static tide and another part reflecting the oscillatory 
properties of the st ar itself. In §13.4.2! we reproduce the 
results presented bv lWillems et a l. (2003) on the orbital 
and spin evolution timescales due to dynamic tides for an 
eccentric binary hosting a 5 Mq MS star and a neutron 
star. We conclude in §|i| CAFein relies on the GNU Sci- 
entific Library (GSL) both for handling the operations 
with matrices and for the integration of the stellar pul- 

^ SPB = slowly pulsating B-type stars 
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sation equations described below. 



yields 



2. COMPUTING NON-ADIABATIC STELLAR 
PULSATIONS WITH CAFEIN 

Before describing in detail CAFein, we give a brief sum- 
mary of the equations governing non-adiaba tic and non- 
radial stell ar pulsations. We refer , e.g., to lUnno et al.l 
(|1989l ) and lGautschv fc Said (|1995n for a detailed deriva- 
tion. 



2.1. The Equations Governing Non-Adiabatic and 
Non-Radial Stellar Pulsations 

The equations governing the non-adiabatic and non- 
radial stellar oscillations are the equations of mass, mo- 
mentum, and energy conservation 



f + V.(pu) = 

d 



p[ — + u ■ V ) u = -Vp - /9V$ 



pT ( — -f u ■ V ) 5* = PEN - V • Ffl. 



(1) 
(2) 

(3) 



To these equations one must add the equations of Pois- 
son and radiative diffusion 



V^* = 4nGp 
„ 4ac 



3ftp 



T^VT. 



(4) 
(5) 



Here p is the density, p the pressure, T the temperature, 
u the fluid velocity, S the specific entropy, $ the gravi- 
tational potential, eN the nuclear energy generation rate, 
Ffi the radiative energy flux, a the radiation constant, 
c* the speed of light, and k the opacity. For simplic- 
ity, electro-magnetic and external forces, viscosity, and 
the coupling of convection and pulsation are neglected 
together with rotation. In the star's frame we take the 
equilibrium state to be spherically symmetric and assume 
that the spatial and temporal part of the small pertur- 
bations can be written in Eulerian form as f'{r,t) = 
/'(r)Y^™(6', (^)e*'^* and similarly for the Lagrangian per- 
turbations, denoted by 6f. Here a is an eigenfrequency 
(below we will denote with oj = \J a-^R^{GM)~^ its di- 
mensionless counterpart), Y^{9,(l)) is a spherical har- 
monic, / is the harmonic degree, and m the azimuthal 
order. 

Following the standard procedure, we then apply a 
small perturbation to the unperturbed star. Perturbing 
and linearizing the basic Eqs. (U)-® and introducing 
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(14) 



Eqs. 



(|9|)- (fT4|) represent the system of equations de- 
scribing non-adiabatic and non-radial stellar pulsations. 
Here, ^ is the displacement of a fluid element from the 
unperturbed position (^r and are its radial and or- 
thogonal components, respectively), g is the local gravity, 
oJ^ = a^R^{GM)~^ is the dimensionless squared eigen- 
frequency, the radiative luminosity, and C4 the ratio 
of thermal to dynamical timescale (rth/Tdyn)- We will 
see below that the latter is used to determine the degree 
of adiabaticity. The remaining terms are summarized in 
Table [T] Recall that the number of radial nodes in 
determines the radial order n of each mode. 

The homogeneous system of Eqs. (|9l)- ([T4| with the 
proper BCs constitute a well-posed eigenvalue problem 
with complex eigenvalue uj. The real and imaginary part 
of the eigenvalue lu (wr and wi, respectively) represent 
the oscillation frequency and the linear growth (wi < 0) 
or damping (lui > 0) rate, respectively. 

In the star's interior, the thermal timescale is much 
longer than the oscillation period and the oscillation be- 
haves almost adiabatically. Therefore, one of the inner 
BCs may be chosen by considering that the entropy is 
constant during a single oscillation (^5* = 0). The other 
BCs at the center are given by the equations of Pois- 
son, mass, and momentum conservation requiring that 
{p'/p + <&'); must be regular at the center: 



ly2 

2/1 2 

2/4 - ly3 = 
2/5 = 0. 



(15) 

(16) 
(17) 



The outer BCs are determined by considering that near 
the surface the Lagrangian perturbation of the pressure 
must vanish, by requiring the continuity of $' and its 
first derivative d$'/d7', and by considering that there is 
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Table 1 

Terms entering the equations governing 
non-adiabatic and non-radial stellar pulsations: Mr 
is the mass contained within a radius r, Cp is the 
specific heat at constant pressure. 
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(18) 

(19) 
(20) 



The homogeneous system of Eqs. (P|- (|14p greatly sim- 
plifies in the adiabatic case, since the terms involving the 
perturbation of the entropy (2/5) and radiative luminosity 
(ye) are neglected. Furthermore, the remaining system of 
four Eqs. (Pl)- (|12p forms an eigenvalue problem with real 
eigenvalue w. In this limit, we require the pressure and 
density to vanish at the stellar surface and Eq. (|18l) re- 
duces to 2/1 — ?/2 + 2/3 = 0. Here we note that neither the 
equations nor the BCs involve the azimuthal order m, 
therefore the eigenvalue is (2^-|-l)-fold degenerate with 
respect to m. 

The radial distribution of the modal families inside a 
star is determined by the run of the Briint-Vaisala (N) 
and Lamb {h\) frequencies, as they characterize the local 
vibrational properties of a star. The Lamb frequency is 



the inverse of the horizontal sound-crossing timescale 

1(1 + \)cl 



(21) 



where Cg = \/Tip/ p is the isentropic sound speed. The 
Briint-Vaisala frequency is the frequency of buoyancy os- 
cillations 

1 dlnp dhip\ 



9 



Ti dr 



dr J 



(22) 



In this work we follow the prescription from 
IBrassard et al.l ()1991f ). which accounts for the buoy- 
ancy due to the gradient in composition. The high 



frequency oscillations (w 



> 



^1 ' 



N'^) have locally the 



characteristics of acoustic waves (p— modes). The low 



frequency oscillations (a; 



< 



^1 ' 



N^) have locally the 



characteristics of gravity waves ((/—modes). 

2.2. The Riccati Method 

Here we describe the numerical method we have im- 
plemented in CAFein to solve the system of Eqs. (O- 
ifT^. This so-called Riccati m ethod, a s intr oduced by 
IScottI (|1973D and extended by IDavevI (|1977f) . was ap- 
pli ed for the first time to the stellar pulsation problem 
bv lGautschv fc Glatzell (|1990aD . 

Commonly, the system of differential Eqs. (|9t- (|14l) is 
solved using relaxation and shooting schemes. These nu- 
merical methods have been successfully applied to solve 
the adiabatic pulsation problem, or to investigate stel- 
lar pulsators which behave almost adiabatically. How- 
ever, commonly used relaxation and shooting methods 
suffer major drawbacks for stars w here non-adiabat icity 
becomes significant (e.g. He stars. ISaio et al.iri984[ ). A 
simple shooting method is stable and feasible for sys- 
tems of two differential equations, as it requires the it- 
eration of only one parameter. However, it becomes un- 
manageable for systems of many differential equations 
like the one described here. In the adiabatic case, if 
the Cowling approximation is adopted (i.e. the per- 
turbation of the gravitational potential is neglected) 
the stellar pulsation problem reduces to a system of 
two ordinary differential Eqs. ^ and ([TU)) and shoot- 
ing methods are normally adopted. The solution thus 
obtained provides an initial guess for the more gen- 
eral non-adiabatic stellar pulsation problem, which is 
enerally solved usin g a Henyey-type relaxation method 
Henvev et"aIlll964D to avoid the cumbersome iteration 
of multiple parameters described above. However, such 
an application of a relaxation scheme will be success- 
ful only if the adiabatic solution is sufficiently close to 
the non-adiabatic one. For cases in which the adia- 
batic and non-adiabatic solutions differ significantly (e.g. 
He stars), the relaxation method often does not con- 
verge and fails to detect certain classes of modes. Be- 
cause of the need for a numerical scheme which does 
not rely on any initial guess es of the adiabatic problem, 
IGautschv fc Glatzell (|1990af ) applied the Riccati method 
to the stellar pulsation problem. This method has 
been extensively and succe ssfully applied to a variety of 
stella r and WD pulsators ("Gau tschy fc Glatze]||1990allbl. 



T99lt iGlatzel fc Gautschy 1991; 



IM?; 'Gautschv & Ldfllcr 1996; 



Glatzel fc KiriakidisI 
Gautsc hv et al.l 119961 : 



iSchenkcr fc Gautschv.1998; .Lofllcr..200Q ). 
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In what follows we provide a detailed explanation of the 
Riccat i method, which we have im ple mented closely fol- 
lowindGautschv fc Glatiell (|1990aD and lTakata fc LofHeii 
(|2004f ). The Riccati method differs from commonly used 
shooting and relaxation techniques mainly in the equa- 
tions that have to be solved, as from a technical point 
of view it is really a shooting method. According to the 
Riccati method the linear first-order ordinary differen- 
tial system describing a boundary eigenvalue problem is 
transformed into a numerically stable, non-linear initial 
value problem. This initial-value problem is then solved 
using a shooting method, where the eigenfrequency is the 
only shooting parameter to be iterated. 

In what follows the subscripts "R" and "I" denote the 
real and imaginary part of complex quantities, respec- 
tively. 



2.2.1. The Riccati Equation 

We start by writing the original system of Eqs. 
in the form 



dr 



A B 
C D 



-m 

(23) 



If the number of elements in the vector y is A^, then M 
is a square matrix of size N x N, while A, B, C, and D 
are square matrices of size J, with N = 2 J. Next, we 
introduce two vectors uric and vric of size J which store 
the first and last J components of y, respectively: 



URi, 
VRi, 



Then, equation (j23l) can be rewritten as: 



d_ 

dr 



URic 
VRic 



A B 
C D 



URic 
VRic 



(24) 



(25) 



From Eq. (j25p we can derive two separate equations for 
URic and VRic 

C^URic 



dr 
dr 



= AURic + BVRic 
= CURic + DVRic 



(26) 
(27) 

(28) 

from Eqs. (1^51) and (P7)) it is straightforward to show that 
R satisfies 



Defining the Riccati matrix R as 



URic 



Rv 



Ric 



dR 

dr 



B + AR RD RCR. 



(29) 



Equation ([^5]) denotes the new system of non- 
homogeneous and non-linear differential equations that 
will be integrated instead of the original stellar pulsation 
problem. Equation has to be solved together with J 
homogeneous BCs at both extrema of the integration in- 
terval. The most general form of the BCs can be written 
as 

PuRic = QvRic (30) 

where P and Q are J x J matrices. As either URic or 
VRic can be considered as arbitrary, Eq. (j30l) uniquely 
determines R. Given the form of the non-adiabatic BCs 



at the star's center ([Tn])-(IT71) and surface ([THl)-(Iini), we 

chose the vectors URic and VRic as 



URi, 




VRi, 




(31) 



The non-adiabatic initial conditions on the Riccati ma- 
trix R at the star's center and surface thus become 



Rc 



Rs 




V+zi 



y(i+2Vadzi) 
2(y+zi) 



V+Z2 

V+zi 

-(' + 1) 

-Z2-2VadV(z2-2l) 
2(y+zi) 



(32) 



(33) 



where the subscripts "c" and "s" denote the center and 
the surface, respectively, and for ease of notation we have 
introduced zi — ^(Z-|-l)/w^ — 4— and Z2 — l{l + l)/uj^ — 
l — l. In the adiabatic case, the BCs on the Riccati matrix 
reduce to 



Rc 



Rs 







1 -1 

-(i + i) 



(34) 



(35) 



2.2.2. The Calculation of the Eigenfrequencies 

For a given frequency w, Eq. (|29p is integrated twice: 
a first integration is performed outward from the star's 
center to a conveniently chosen fitting point rat with ini- 
tial conditions ([32]) yielding matrix R°"*(r). A second 
integration is performed inward from the star's surface 
to rfit with initial conditions ([33]) yielding matrix R'"(r). 
The frequency w is an eigenfrequency, when the eigen- 
functions ?/i (i = 1 — > 6) are continuous at the fitting 
point: 



URic 
VRic 



URic°"*(rfit) 
°"'(rfit). 



VRi, 



(36) 
(37) 



Eqs. (|36|) and ((37|) are equivalent to 



[R'^rtit) - R° 

A necessary condition for Eq 
solution yields 



(rfit)]vRic = 0. (38) 
to have a non trivial 



det[R'"(rfit)-R°"*(rfit)] =0. 



(39) 



Since uj is complex, we first scan the parameter space 
in ojR setting uJi = 0. At this stage we find the inter- 
val in wr across which the real part of expression (|39p 
crosses zero. Next, the values of ojr at the extrema 
of this interval are taken together with the values of 
det[R'"(rfit) — R°"*(rfit)] as initial guesses for the iter- 
ation of the exact eigenfrequencies. During this itera- 
tion we use a complex secant method to find the ex- 
act values of wr and uJi for which both the real and 
im aginary part of conditio n (|39l) are satisfied. Follow- 
ing l^^ta^XHffleB ()2004[ ) we chose the fitting point rat 
based on the behavior of the Briint-Vaisala and Lamb fre- 
quencies. Specifically, for a given frequency wr we pick 

: u!^. Scanning the 



the innermost point where N'^ or 
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parameter space in w following the procedure described 
above yields all eigenfrequencies. 

2.2.3. Reimbedding 

During the calculation of the eigenfrequencies fii l2.2.2p . 
VRic might vanish and therefore R becomes singular (re- 
call Eq. ([281) '). Ill this case, we avoid the singularity by 
re-defining upic and vric according to the so-called "re- 
imbedding" procedure, as described by ,Takata & Loffie^ 



(|200l . The new u'^^^^ and 



'Ric 



are a transformation of 



the original URic and vric and are given by 



''Ric 



URic 
VRic 



Too Toi 
Tio Til 



URic 
VRic 



(40) 



where T is a square N x N matrix, while Tij are square 
J X J submatrices (recall N = 2 J). For the transformed 
Riccati matrix R' we can write 

R' = (TooR + Toi)(TioR + Tii)-i, (41) 

choosing matrix T such that 

det(TioR + Tii) ^0. (42) 

The transformed Riccati matrix R' is still a solution for 
Eq. (pn)) . provided that A,B, C, and D are replaced by 
their prime (') counterparts. 



A' B' 
a D' 



T 



A B 
C D 



(43) 



We apply this procedure whenever the Eucledian norm 
of R (||R|| = 



jj=u'-i3' g06S above a certain 
pre-defined value, switching to the permuted R' with 
the lowest ||R'||. Here we note that for the fully non- 
adiabatic case, the integration variables (i = 1 — )> 6) 
are complex. In this configuration, URic and vric have 
size 6, and we chose them so that 



URic = 



/ 2/13 \ 

2/4,R 
2/5,R 



2/4.1 

\y. 



URic,R 
URic. I 



'5.1 



VRic = 



/ 2/2,R \ 

2/3,R 

2/6,R 

2/2J 
. 2/3J , 
\ 2/6J / 



VRic,R 
VRic, I 



(44) 



The resulting R has a size 6x6 and the search for the 
permutation yielding R' with the minimum norm would 
result in long computational times. However, given that 



URic,R 

URic, I 



Ri R2 

R3 R-4 



VRic,R 
VRic, I 



(45) 



we minimize the numbers of trials finalized to find the 
permutation yielding R' with the minimum norm by only 
permuting the (dominant) real part of R (Ri and R4) 
and applying that same permutation to the imaginary 
part (R2 and R3). In the adiabatic case, instead, only 
Eqs. (|9))-p2|) are solved (the terms containing 2/5 and tjq 



Table 2 

Units of Physical Quantities. The total mass, 
radius, and surface luminosity are denoted 
with M, R, and L, respectively. 



Quantity 



Unit 



Radius (r) 
Mass (m) 
Pressure (P) 

Radiative Luminosity (I/r) 

Frequencies {L'^,N'^,uP) 

Density (p) 

Gravity (g) 

Temperature (T) 

Sound speed (cg) 

Specific heat, entropy (cp, S) 



R 

M 

GM'^/{4ttR^) 
L 

GM/R^ 
M/{'inR'^) 
GMjR^ 
GMIi R^^^R) 
^GM/R 

Rs^as 



are neglected) and the eigenfrequency is purely real. This 
reduces the size of the Riccati matrix to 2 x 2 and we scan 
on all possible permutations to minimize ||R'||. 

2.2.4. The Calculation of the Eigenfunctions 

Once the eigenfrequencies of interest have been deter- 
mined, the calculation of the eigenfunctions for a par- 
ticular eigenfrequency proceeds as follows. Eq. (|29)) is 
integrated as described in § 12.2.21 and the components of 
R (-Rij) are stored together with the permutations ap- 
plied. Once Rij have been evaluated both for R'"(r) 
and R°"*(r), the eigenfunctions are readily calculated by 
solving 

^ = (CR + D)vRie (46) 
ar 



together with Eq. (1^51) . It is straightforward to derive 
Eq. dMl) from Eqs. dS?]) and ([Ml)- Eq. (gS]) is integrated 
twice, from the fitting point rgt to the star's center and 
from Tfit to the surface. At this stage, we do not inte- 
grate R again, but we interpolate the already calculated 
Rij using a third order polynomial (Steffcn 1990) and 
apply the same permutations used during their calcula- 
tion. The initial value of vric is given by a non-trivial 
solution of Eq. ((55)) . Clearly, during the integration from 
Tfit to the center (ret to the surface) R°"* (R™) must 
be used. This numerical scheme where vric and R are 
not integrated together and are integrated in the oppo- 
site d irections is necessary for numerical stability ((Sloanl 
IT971 . 

2.3. Examples of Free Stellar Oscillations and Tests 

In this section, we present a series of tests for 
CAFein both internal and against published results 
in the literature. For the latter, we compute the 
eigenfrequencies of a polytropic model and compare 
them with results in the literature wh i ch rel y both on 
the Riccati method (" Takata fc Loffieri 12004 hereafter 
TL04) and on standard relaxation and shoo ting tech- 
niques (jChristensen-Dalsgaard &: Mullanlll994 hereafter 
CDM94). In what follows, we refer to the eigenfrequency 
of the purely adiabatic problem with WAd, while we de- 
note with cj the non-adiabatic eigenfrequency. The sub- 
scripts "R" and "I" have their usual meaning. 

2.3.1. Testing CAFein 's numerical accuracy on a 
Polytrope: the Adiabatic Case 
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Figure 1. Test on CAFein's performances if the required inte- 
gration accuracy (both relative and absolute), ||R||iimiti and res- 
olution adopted during the scan of the parameter space in uj are 
varied. For the latter, we consider a randomly chosen interval 
[a;Ad-0.1, WAd+0.1] across each eigenfrequency. Here we set I = 2, 
hut the same test performed on the I = 3 eigenfrequencies yielded 
similar results. The relative difference is the absolute value of the 
difference between our fiducial eigenfrequencies and the ones calcu- 
lated by changing the parameters mentioned above, divided by our 
fiducial values. Top: the integrator accuracy and ||R-||iimit are 
kept fixed, while the eigenfrequency scan resolution (# steps) is 
varied; Middle: the integrator accuracy and eigenfrequency scan 
resolution are kept fixed, while ||R.||ii,nit is varied; Bottom: same 
as the middle panel, but for a lower integrator accuracy (note that 
an accuracy of 10~^^ and ||R||iimit = 10 determine our fiducial 
values). Some of the data overlap and are not visible. 

In this section, we apply CAFein to a polytropic model 
to compute its eigenfrequencies. We then test the nu- 
merical accuracy of our calculation by changing some of 
the relevant parameters entering the computation of the 
Riccati matrix and analyzing how the eigenfrequencies 
vary. In what follows, when we mention the integrator 
accuracy adopted in our calculation, we refer both to the 
absolute and relative accuracies. 

We create a polytrope with index n — assuming 
an ideal gas with Fi ~ 5/3. Below, we indicate the 
polytropic index with Up to avoid confusion with the in- 
dex denoting the radial order of a mode. We first solve 
the Lane-Emden equation using a variable step 4th or- 
der Runge-Kutta integrator with an accuracy require- 
ment of 10~^^. We then iterate the integration until we 
reach a resolution between two consecutive mesh points 
of Ar < 5 X 10^^ (here we are following TL04, but adopt- 



ing a slightly higher resolution). We conveniently pass 
onto dimensionless quantities by expressing the physi- 
cal quantities in the units listed in Table [21 This ta- 
ble also lists the variables that will be used in the non- 
adiabatic regime. During the calculation of the adiabatic 
eigenfrequencies, we use again a variable step 4th order 
Runge-Kutta integrator and we interpolate the various 
polytropic pa rameters using a third order polynomial 
(jSteffenl 119901 ). As an example, we report some of the 
calculated eigenfrequencies uj Ad in Table [3l Here we set 
the required integration accuracy to 10^^^, the limit on 
the Eucledian norm of R (||R|| limit) used for rcimbcd- 
ding (ii l2.2.3p to 10, and perfomed 1000 integrations in w 
on a randomly chosen interval [wAd-O.l, WAd+O-l] across 
each eigenfrequency. In what follows, we will refer to the 
eigenfrequencies listed in Table [3] as our "fiducial" ones. 

Next, we test the stability of CAFein with respect to 
the parameters mentioned above. We re-compute the 
same eigenfrequencies for different integrator's accura- 
cies (10-'' and 10-12), ||R||iimit (10, 100, and 1000), and 
number of integrations across each eigenfrequency (50, 
200, and 1000), fixing the width of each integration in- 
terval as described above, and compare them with our 
fiducial values. The outcome of this test is summarized 
in Fig.[l] where it is clear that our results are sensitively 
affected only by the resolution adopted during the scan 
of the parameter space in uj (i.e. the eigenfrequency scan 
resolution, see top plot). Varying ||R||iiinit while keeping 
the number of iterations and integrator accuracy fixed, 
the relative difference with our fiducial eigenfrequencies 
is between IQ-^'^ — 10^^^ (middle plot in Fig.[T]). Decreas- 
ing the integrator accuracy by 3 order of magnitudes, 
the relative difference with our fiducial eigenfrequencies 
is < IQ-* (bottom plot in Fig. [1]). The same test per- 
formed on the I = 3 eigenfrequencies yielded similar re- 
sults. 

Once the eigenfrequencies are known, the eigenfunc- 
tions can be readily calculated as described in § 12.2.41 
The e igenfunctions should be orthogonal. In fact, as 
[Fuller fc Lai (2011) point out, the numerical determina- 
tion of an eigenfunction might be contaminated by other 
eigenfunctions 

(la)„um = hai^ + ho^o + hiii + ■■■ (47) 

where the subscript a denotes the order of the mode and 
the displacement from the equilibrium position for the 
mode a can be written as 

€a = Kr,a(r)e, + a,c. (r)eh V] (48) 

The terms with subscript "0" in Eq. (|47p refer to the 
/-mode and the various coefficients are given by 

/Ji=<€i|$„ >= / r2p[er,ier,a+K^ + lKh,i^h,Jdr (49) 
Jo 

where the various eigenfunctions in Eq. (j49|) are the ones 
calculated numerically. To derive Eq. (|39|). we used the 
normalizati o n of the spherical harmonics as given by 
lUnno et all (fl989l ): 

/ / Yr{e,cjj)Yir-' (9,^) sine d9dcj)^Swd,nm' (50) 
JQ Jo 

where 6w and Smm' are the Kronecker deltas. Since the 
/—mode gives the dominant contribution in Eq. (|47p . we 
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Table 3 

Eigenfrequencies for a polytropic model with np = 3. Here 
we used an accuracy requirement for the integrator of 10~^^, 
l|R-i!limit = lOj ^iid 1000 integrations in w on a randomly 
chosen interval [ujj^^-O.l, uij^^+O.l] across each 
eigenfrequency. See Fig. [T] for a test on the numerical 
accuracy of the eigenfrequencies calculated in this work and 
Fig. [2] for the relative differences between our results and the 
ones presented by TL04 and CDM94. For TL04 the values are 
taken from Table 1 of their paper, while for CDM94 we used 
Table 2 and 4 of their paper. 



Ref. 


mode 






This Work 
TL04 
CDM94 


P40 


2.7777508707x103 
2.7777508750x103 
2.7777509770x103 


2.8357237896x103 
2.8357239610x103 
2.8357238770x103 


This Work 
TL04 
CDM94 


P30 


1.6220938771x103 
1.6220939110x103 
1.6220938720x103 


1.6652407610x103 
1.6652408790x103 
1.6652410890x103 


This Work 
TL04 
CDM94 


P2() 


7.7357674655x10^ 
7.7357675130x10^ 
7.7357672120x10^ 


8.0228405811x10^ 
8.0228406290x102 
8.0228405760 xlO^ 


This Work 
TL04 
CDM94 


PlO 


2.3362818835x10^ 
2.3362820270x10^ 
2.3362818910x10^ 


2.4860043045x10^ 
2.4860044040 xlO^ 
2.4860043330 xlO^ 


This Work 
TL04 
CDM94 


Pi 


1.5263662310x10^ 
1.5263662338x101 
1.5263660431x101 


1.8443609723x10^ 
1.8443608440x101 
1.8443605420x101 


This Work 
TL04 


/ 


8.1753397221 
8.1753397230 


9.4137919393 
9.4137926170 


This Work 
TL04 
CDM94 


91 


4.9145734152 
4.9145734160 
4.9145731920 


6.7669725650 
6.7669720220 
6.7669711110 


This Work 
TL04 
CDM94 


gio 


3.2249531558x10^ 
3.2249534130x10^ 
3.2249531150x10^ 


5.8751831508x101 
5.8751833439x101 
5.8751821518x101 


This Work 
TL04 
CDM94 


920 


9.7498882680E-02 
9.7498878837E-02 
9.7498863935E-02 


1.8574303847x101 
1.8574304879x101 
1.8574303389x101 


This Work 
TL04 
CDM94 


930 


4.6535316784x10^ 
4.6535320580x10^ 
4.6535316855x10^ 


9.0075142564x10^ 
9.0075135231x102 
9.0075127780x10^ 


This Work 
TL04 
CDM94 


940 


2.7186954918x10^ 
2.7186956257x10^ 
2.7186954394x10^ 


5.3050607114x10^ 
5.3050611168x102 
5.3050607443 xl02 



normalize the eigenfunctions so that ha =< ^^\^^ > — 1 
and we take {$,a)num to accurately represent the actual 
^a if 1^0 1 *C 1- The values of the coefficient ho for the 
/—mode and the first five p— and modes of harmonic 
degree I = 2 for the polytropic model considered are re- 
ported in Table m The results show that orthogonality 
is satisfied. 

Recall that CAFein has been developed to investigate 
dynamic tides in close binaries. Our focus on the har- 
monic degree I = 2 will become clear in § 13. 1[ where we 
introduce the tide-generating potential. 

2.3.2. Comparing our Fiducial Polytrope Eigenfrequencies 
to Polytrope Results in the Literature 

In this section, we compare our fiducial eigenfrequen- 
cies calculated in the previous section with results in 
the literature which rely both on the Riccati method 
(TL04) and on standard relaxation and shooting tech- 
niques (CDM94). Our fiducial eigenfrequencies are listed 
in Table |31 together with the eigenfrequencies calculated 
by TL04 and CDM94. In Fig. [2] we show the rela- 
tive difference between our fiducial values and TL04 and 
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Figure 2. Relative difference between our fiducial eigenfrequen- 
cies and the one presented by TL04 and CDM94. As in TablefS] we 
consider a polytropic model with np = 3. Top: 1 = 2; Bottom: 
I = 3. The relative difference is the absolute value of the differ- 
ence between our fiducial eigenfrequencies and the ones reported 
by TL04 or CDM94, divided by the fiducial values. For compari- 
son, we denote with "x" the relative difference between TL04 and 
CDM94. 



Table 4 

Orthogonality of the I = 2 eigenfunctions 
for a Tip = 3 polytropic model. The 
coefficient \ho\ for the /—mode is 1. 



mode 








mode 




l^ol 




P5 


1.4 


X 10" 


6 


91 


5.4 


X 10" 


6 


Pi 


8.4 


X 10" 


-7 


92 


4.3 


X 10" 


6 


P3 


2.3 


X 10" 


-6 


93 


2.3 


X 10" 


6 


P2 


5.5 


X 10" 


-6 


94 


1.3 


X 10" 


6 


Pi 


3.2 


X 10" 


-5 


95 


1.3 


X 10" 


6 



CDM94 results (with filled circles), and the relative dif- 
ference between TL04 and CDM94 (with "x"). The 
relative difference between our fiducial values and TL04 
(CDM94) is between - 10"" - IQ-^ (- 10"* - lO"^) 
both for I = 2 and ^ = 3. The upper end of this in- 
tervals agrees with the relative difference between TL04 
and CDM94. In particular, the better agreement with 
TL04, results in a relative difference between our fidu- 
cial values and CDM94 that is nearly the same as the one 
between TL04 and CDM94 (orange filled circles overlap 
with "x"). This behavior occurs for both the I — 2 and 
I — 3 eigenfrequencies. Among the factors that might 
be contributing to the small differences between the var- 
ious results presented here are: the different accuracies 
adopted during the calculation of the polytropic model 
among the different studies, the accuracy adopted dur- 
ing the calculation of the eigenfrequencies, and round-off 
errors, as suggested by TL04 ( § 5 of their paper). 
Hence forward we focus on the harmonic degree 1 = 2. 
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2.3.3. Testing CAFein on a Zero Age Main Sequence Star: 
the Non-Adiabatic Case 

In this section we test the behavior of CAFein in 
the non-adiabatic regime by applying it to a Zero Age 
Main Sequence (ZAMS) model in the (3 Cephei/SPB 
region of the HR diagram. Our goal here is only to 
verify the reliability of CAFein in identifying unsta- 
ble modes; a detailed study of (3 Cepheis and SPBs 
is beyond the scope of this work. We refer to, e.g., 
IGautschv fc Sa io (1995, 199|) for a review of pulsating 
stars and to Cox et al. (1992), Moskalik & Dziem bowski 
( 1992j) . Kiriakidis et al.. (1992,,). and Dziemb owski et al. 
()1993l ) for detailed investigations targeting (3 Cephei and 
SPB variables. 

Insofar as the stellar mod e l ado pted here is con- 
cerned, we follow ISaiolTCo^ (Il980l) (hereafter SC80), 
who investigated /3 Cepheis near the MS and originally 
found all the models investigated to be stable. Stabil- 
ity was due to th e use of opacity formula e prior to the 
one proposed by iRogers fc Iglesiad (I1992D. Only after 
the new OPA L opacitv tables ("Rogers fc lglesiasi 119921: 
Iglesias et al.l [19 92: Seaton et al. 1994; Kiriakidi s et al.l 
1992t iMoskalik fc Dziembowskil 119921) were introduced 
did it became clear that the excitation mechanism re- 
sponsible for pulsations in these stars is the so-called 
K-mechanism due to an opacity bump in the heavy ele- 
ments. This is in contrast to stars located in the classical 
instability strip whose oscillations are driven by the k- 
mechanism due to partial ionization of H and He I and/or 
He II. 

Following SC80, we use MESA (jPaxton et al.l[20TTI ) to 
create a ZAMS model of 7 Mq at metallicity Z = 0.03 
and X = 0.7. We then increase the number of mesh 
points by i nterpolating the model with a third order 
polynomial ()Steffenlll990D . to reach a resolution between 
two adjacent mesh points of Ar < 5 x 10~^. Our 
stellar model has a luminosity, effective temperature, 
and radius of log(L/LQ) = 3.249, log(Teff/K) = 4.309, 
and log(i?/Ro) = 0.530, respectively. The H abun- 
dance at the center is 0.67. For comparison, the same 
properties for one of the models used by SC80 are 
log(i/Lo) = 3.246, log(Teff/K) = 4.311, and log(i?/RQ) 
= 0.523, while the H abundance at the center is 0.7. 
We conveniently pass onto dimensionless quantities by 
expressing the physical quantities in the units listed in 
Table [H From these quantities we then derive the param- 
eters introduced in § 12.11 and the one listed in Table [TJ 

To have a sense of where the modes can propagate in- 
side the star and whether non-adiabatic effects are signif- 
icant, we calculate a few adiabatic eigenfrequencies and 
place them on the propagation diagram (see Fig. [J]). This 
diagram shows that modes can propagate all the way 
to the surface, while 5— modes are more confined towards 
the star's interior, especially the low-order modes. To in- 
vestigate the importance of non-adiabatic effects for this 
stellar model and for the modes of interest, we consider 
the ratio of the star's thermal to dynamical timescale 
(^th/Tdyn)- Siucc the fundamental oscillation timescale 
(defined as the travel time of a sound wave from the cen- 
ter to the surface) is of the same order of magnitude as 
Tdyn, non-adiab atic effects are re levant when Tth/Tdyn is 
small (see e.g. lUnno et al.l 119891) . The dashed-dot line 
in Fig. |3] shows that even though Tth >> Tdyn through 
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Figure 3. Propagation diagram and degree of non-adiabaticity 
for the 7 Mq ZAMS star described in the text. Left y-axis: the 
solid lines denote the Briint-Vaisala (A'^) and Lamb (Lj) frequencies 
squared, while the dots represent the zeros of the radial part of the 
eigenfunctions of the modes (from the top to the bottom) p2,pi,gi- 
gg. The harmonic degree was set to i = 2 and the units in Table [2] 
are used. Right y-axis: the dashed-dot line shows ratio of the 
thermal timescale (rtjj) to the dynamical timescale (r^yn). Non- 
adiabatic effects become significant when the two timescales are 
comparable (see text). 



Table 5 

Non-adiabatic I = 2 eigenfrequencies for 
the 7Mq ZAMS described in the text. 



mode 








UlJ 






P2 


5.16434 




1.01 


X 


10- 


-4 


Pi 


4.01821 




-6.17 


X 


10- 


-7 


f 


3.18180 




-2.14 


X 


10" 


-7 


gi 


1.86538 




-1.71 


X 


10" 


-7 


92 


1.27920 




-3.75 


X 


10" 


-7 


93 


9.67362 X 10- 


1 


-5.84 


X 


10" 


-7 


9i 


7.71946 X 10- 


1 


-7.91 


X 


10" 


-7 


95 


6.39703 X 10- 


1 


-9.65 


X 


10" 


-7 


96 


5.45513 X 10- 


1 


-1.06 


X 


10" 


-6 


97 


4.76716 X 10- 


1 


-9.75 


X 


10" 


-7 


9S 


4.25366 X 10" 


1 


-4.06 


X 


10" 


-7 


99 


3.85168 X 10- 


1 


1.14 


X 


10- 


-6 



the bulk of the star, non-adiabaticity becomes relevant 
approaching the surface, where p-modes can propagate. 

Next, we follow the procedure outlined in § 12.2.21 and 
calculate the non-adiabatic eigenfrequencies. As the in- 
clusion of non-adiabatic effects renders the stellar pulsa- 
tion equations stiff, at this stage of the calculation we use 
a variable step implicit Bulirsch-Stoer integrator with an 
accuracy requirement of 10~^^. The results are summa- 
rized in Table [5l From a test on CAFein's numerical sta- 
bility like the one performed in § I2.3.1[ we find that the 
calculated eigenfrequencies are accurate as long as the re- 
quired integrator accuracy (both absolute and relative) 
is < 10^^^. As a negative uji denotes an unstable mode, 
we can see the excitation of the modes which lie in the 
transition region bet ween q— and p— rnodes. as expected 
for 13 Cepheis (e.g. IGautschv fc Salol 119961 ) and SPBs. 
The latter are generally understood as an extension of 
the (3 Cephei instability towards longer periods (smaller 
frequencies), as their observed pulsation periods are due 
to the excitation of g— modes. As a reference, the peri- 
ods for the pi— and gg— modes listed in Table [5] in days 
are ~ 0.07 and ~ 0.6, respectively, which is consistent 
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with the range of oscillation periods observed for these 
kind of stars. Here we also note that the magnitude 
of LUi (and therefore non-adiabaticity) is negligible for 
(7— modes, while it increases by about two order of mag- 
nitudes for the mode. This was expected given the 
trend of Tth/Tdyn shown in Fig. |3l Since the p2— mode 
is the only mode considered here for which dissipation 
is significant, we calculate the non-adiabatic eigenfunc- 
tions for this mode. Recall that the calculation of the 
eigenfunctions is performed in two steps. A first inte- 
gration yields the components of the Riccati matrix and 
the permutations applied, while a second integration uses 
these for the calculation of the eigenfunctions. During 
the first integration CAFein uses a Bulirsch-Stoer integra- 
tor with an accuracy requirement of 10~^^ and it then 
refines the number of mesh points in close to each 
permutation point using a Runge-Kutta integrator with 
the same accuracy requirement. During the calculation 
of the eigenfunctions we integrate Eq. (j^S]) using again 
a Bulirsch-Stoer integrator and the same accuracy re- 
quirement. Following the representation of the /—mode 
eigenfunctions in Fig. 1 of SC80, we show in Fig.|4]some 
of the non-adiabatic eigenfunctions for the p2-inode. The 
purely adiabatic radial component of the displacement 
from the equilibrium position (dashed line) is also shown 
for comparison. As expected, (^r/'')R and its adiabatic 
counterpart are very similar and the entropy perturba- 
tion [both {6S/cp)-[i and {SS/cp)i] increases rapidly to- 
wards the star's surface, where non-adiabaticity becomes 
significant. Similarly to the results presented by SC80, 
the real part of the entropy perturbation presents two 
minima which are located at the peaks of the opacity k. 
For the case of SC80, the peaks are in the He ioniza- 
tion zone. In our MESA model, the lower-temperature 
bump in the opacity is due to the He ionization, while 
the one at a higher temperature is due to photon ab- 
sorption by the L-shell of F e and photoionization f rom 
the K -shell of C, O, and Ne (IRoeers fc IglesiaslfTOOl) . It 
is the fc— mechanism associated with this second bump at 
a temperature of ~ 2 x 10^ K which drives the pulsations 
observed in /3 Cephei variables. 

3. INVESTIGATING DYNAMIC TIDES WITH 
CAFEIN 

Before describing the extension of the Riccati method 
developed to investigate tidally excited stellar oscilla- 
tions, we briefiy outline the basic assumptions adopted in 
this work and introduce the various parameters entering 
the dynamic tides theoretical framework. 

3.1. The Tide- Generating Potential 

We consider a close binary system of stars with masses 
All (primary) and M2 (secondary) orbiting around one 
another in a Keplerian orbit. We assume that the pri- 
mary has a radius Ri and that it rotates uniformly 
around an axis orthogonal to the orbital plane with an- 
gular velocity fli in the sense of the orbital motion, 
while we treat the companion as a point mass. Fur- 
thermore, we assume fii to be small enough so that 
the Coriolis force and the centrifugal force can be ne- 
glected. Under these assumptions, the tides raised by the 
companion can be treated as small forced perturbations 
applied on a spherically symmetric star in hydrostatic 
equilibrium. Following the general procedure, we can 
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Figure 4. Non-adiabatic I = 2 p2-mode for the 7 M© ZAMS 
star described in the text. Some of the eigenfunctions are 
shown as a function of the normalized stellar radius close to 
the surface. The eigenfunctions have been normalized so that 
S,i/t = (Cr/'')R -I- i{£,r/r)i = 1 at the star's surface. Recalling 
that yi = ^v/f and expressions l|6ll-(|8ll, this normalization yields 
for each eigenfunction y, {j/i,ri/r + yi,iyi) / (vl^^ -I- ?;f j) Vn and 
(s/i,rW - !/i,i?/r)/(?/?,r + yl,i) J/I> where the real and imag- 
inary part of yi are evaluated at the surface. Solid and dotted 
lines represent the real and imaginary part, respectively. The adi- 
abatic radial component of the displacement from the equilibrium 
position is also shown for comparison (dashed line). 

express the tide-generating potential in spherical coordi- 
nates r = (r, 9, 4>) with respect to an orthogonal frame 
corotati ng with the star and exp and it in Fourier series 
as (e.g. iPolfliet fc Smeverslll99dt hereafter PS90) 

e^W{v,t) = -e^Yl E E Ci.^.J;^ Yr{e,<f>) 

X exp[j(a„i,fct - kQovbr)] (51) 

where the polar angle 9 is measured from the rota- 
tional angular velocity vector, while the azimuthal angle 
(f) is measured in the orbital plane and in the sense of 
the orbital motion. At time t = 0, the angle cj) = 
marks the position of the periastron of the binary orbit. 
The tide-generating potential is a solution to Laplace's 
equation. The indices l,m, and k in Eq. ()51|) are the 
harmonic degree, the azimuthal number, and the Fourier 
index, respectively. The dimensionless parameter ex = 
{Ri / a)^ {M2 / Ml) measures the ratio of the tidal force to 
gravity at the star's equator, a is the semi-major axis, 
o'm,fe — fcr^orb + m^i is a forcing angular frequency with 
respect to the corotating frame, riorb = 27r/Porb the 
mean motion, r a time at periastron passage, and ci^rn,k 
are Fourier coefficients defined as 

_ (^-|m|)! /fiiV-^ 1 1 

'■"'""(^-flml)! ' ^'\a) (l-e2)'-i/2^ 

X / (1 4- e cosi/)'"^cos(fcM + mu)dv. (52) 
Jo 

Here, P[^{cos9) are Legendre polynomials of the first 
kind, v is the true anomaly and M = f2orb(i — t) the 
mean anomaly. The main properties of the Fou rier co- 
effic ients were desc r ibed b y lSmevers et al.l ()1998[ ). PS90, 
and iWillems et al.l (|2010[ ). Briefly, ci^rn,k are symmet- 
ric with respect to m and k {ci^m,k = ci^-m,-k) and are 
equal to zero for odd values of /-|- |m| since P™(0) = for 
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odd values of 1+ \m\. Furthermore, the binomial theorem 
implies that ci^m.o = 0. For a given orbital eccentric- 
ity, the absolute value of ci^rn,k decreases with increasing 
k, th ough the decrease is slower for high l y eccentric or- 
bits (IWillemsl[200a IWiTlems et all [2003HSmevers"eraII 
119981 ). This implies that the number of ci^rn,k terms with 
non-trivial contributions to the tide-generating potential 
is finite, though it increases with increasing eccentricity. 
Given the dependence of Q_m,fc on (i?i/a)'~^, investiga- 
tions on dynamic tides are often restricted to the terms 
belonging to Z = 2, as they are dominant. 

It is clear from the expansion ((5T|) of the tide- 
generating potential that the tidal action from the com- 
panion induces in the primary an infinite number of forc- 
ing angular frequencies (Tm,fc. The terms associated with 
o'm,k = (the time-independent terms in exW^) give rise 
to static tides, while the terms associated with am,k 
(the time-dependent terms in stW) give rise to dynamic 
tides. 

3.2. The Equations For Tidally Excited Stellar 
Pulsations and The Secular Evolution of the 
Orbital Elements 

If the tides raised by the companion are treated as 
a small perturbation applied on a spherically symmet- 
ric star in hydrostatic equilibrium, the equations de- 
scribing forced stellar oscillations are still derived from 
Eqs. (P)-®, provided that the erW is ad ded to Eg. (|2t 
of m omentum conservation (e.g. PS90, iWillems et al.l 
120101 ). Following the standard procedure, we take the 
unperturbed solution to be axisymmetric and assume 
that the spatial and temporal part of a small perturba- 
tion can be written in Eulerian form as /.^(r, 0, t) — 

Et=2EL^-iEZ-^fUA^)yr(s,^y^-;'' or simi- 

larly for the Lagrangian form, denoted with S. Since 
the tide-generating potential is a solution to Laplace's 
equations, perturbing and linearizing the basic Eqs. ([T])- 
([5|) with the new equation for momentum conservation 
yields, for each set of {l,m,k) in the expansion of the 
tide-generating potential, a system of equations which 
is formally identical to Eqs. (Pt- P^ . with the follow- 
ing modifications. The perturbation of the star's grav- 
itational potential and the tide-generating potential are 
grouped into the total perturbation of the gravitational 
potential defi ned as ^ + t^W (e.g. lZahnlll975L 

PS90, Willem s~eraIll2OT0l) . where (denoted in S [Q 
as $') is the perturbation of the star's potential of self- 
gravitation due to the tidal action of the companion. Fur- 
thermore, the new integration variables retain the same 
form, provided that $' is substituted with ^ . A final 
modification to Eqs. (|9|)- ([T4|) concerns the BCs at the 
star's surface. As the gravitational potential and its first 
derivative must be continuous at r = BC (fTH)) be- 
comes (e.g. PS90) 



y4 + (^ + 1)2/3 + 2/1 eT(2Z ^ l)c;,„,,fc = 0. 



(53) 



Therefore, the introduction of keeps the tidally ex- 
cited stellar pulsation Eqs. (|9|)-(fT4|) homogeneous, but it 
renders the BCs non-homogeneous. Because of the non- 
homogeneous term in Eq. ([S^]) . the solutions to Eqs. 

are proportional to eTCi,m,/c- Furthermore, even 
though the system of equations is complex, the dimen- 
sionless tidal forcing frequency Wm.fc, is purely real (recall 



that J, — f.R^{GM) ^). In what follows, we refer 
to the solution of Eqs. ©-([HI) with BCs ([T5l)-(IT7l) at the 
star's center, and BCs (fTSl) . ([53| . and (f2Ul) at the star's 
surface with "tidal eigenfunctions" . 

From the tidal eigenfunctions, the timescales for the 
secular evolution of the orbital elements and stellar spin 
due to dynamic tides can be readily calculated. 

The evolution of the orbital separation and eccentric- 
ity is due to the primary's tidal deformation, which in 
turn perturbs the external gravitational field and there- 
fore the Keplerian motion of the binary components. En- 
ergy dissipation in the surface layers causes a phase shift 
between the perturbation of the gravitational potential 
and the companion's position. This phase shift results 
in a torque exerted from the secondary on the tidally de- 
formed primary, which affects the primary's spin. The 
rates of secular evolution for a, e, and fti are given by 
(e.g. IWillems et aTllMol) 



! = 1 m=-l fe = ^ ^ 

,m,k feG ™ Ae) 

4 I oo , „ \ 1+3 



Stt Ma 



forb 



a 

'(3) 



1 = 1 m = -l ft=0 
\Fi ,m,k l^i^'y I, m,kG I fc ('^) 

>orb \Mr+M2) Ml h 



(54) 



(55) 



1 = 1 m = -l k=0 



X sin 7i.m,feG,*ti,fc(e) 



(56) 



where Ii is the star's moment of inertia and Eq. ([56| is 
derived assuming solid-body rotation and that the tidal 
deformation does not affect Ji. In the above equations, 
the dimensionless Fi^m.k measure the response of the star 
to the various tidal forcing frequencies and are given by 



F, 



Rl '^Lm,kiRl) 
GMi eTQ,m,fc 



1 



(57) 



These coefficients because 
'^i,7n,k{Ri) erci^m.k 13. 2p . In the units of Table [H 
expression (|57|) reduces to 



F, 



(2/3)i,m,fc(l)g(l) 



1 



^\F, 



l^m^k I 



(58) 

where the last equality comes from the complex na- 
ture of the tidal eigenfunctions. For the various prop- 
erties of symmetry obeyed by \Fi^m,k\ and for the def- 
inition of Ki^m.k (not to be confused with the opacity). 



'^lm.fc(e). fi\m.fc(e): and G);^ fc(e), a nd their properties 
we refer to IWillems et al.l (|2003ll2010D . Here we just note 
that GJ ^ ^(e) are all zero for a binary with a circular or- 
bit. Eqs. ([M)) - ([5^ take the same form as the equations 
for the rate of secular change of orbi t al separatio n ec- 
centricity, and spin derived bv ' ZahnI (|1977l . I1978D . iHutI 
(|1981i) , and iRuvmaekera (,1992^ , in the limiting case of 



.(3) 



(4) 
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weak damping and small forcing an gular frequencies (see 
Appendix D of lWillems et al.ll201Cll for a derivation). 

In what follows, we omit the subscripts /, to, and k 
from the components of the tidal displacement field and 
the perturbed stellar structure quantities and we denote 
the tidal forcing frequency by wt- 

3.3. Extending the Riccati Method to Forced Stellar 
Oscillations 

As described in § 12. 2[ the Riccati method relies in go- 
ing from a homogeneous system of ordinary differential 
equations to a non-omogeneous one. However, as ex- 
plained in § 13. 2| even though the equations describing 
the tidally excited stellar oscillations are homogeneous, 
the same is not true for the BCs at the star's surface. 
We make the Riccati method viable for investigating dy- 
namic tides by introducing two new variables yj and ys, 
such that BC (|53l) at the star's surface becomes homo- 
geneous 

y4 + {l + l)y3 + yi + eT(2/ + l)ci„n,ky8 = 0. (59) 

The introduction of two variables instead of one is re- 
quired to keep all the matrices entering the Riccati 
method square. We take yr and ys to be a solution of 
the following differential equation 



dyv dya 
r—r- — r- 







(60) 



dr dr 

with BCs at the center and at the surface given by 

yr = 2/8- (61) 

Once the tidal eigenfunctions are determined, we nor- 
malize them so that y% — 1 at the star's surface. This 
choice of normalization causes Eq. ([59l) to reduce to the 
original BC ((53|). Here we note that changing the form 
of Eq. (|60p does not affect our results, but it can affect 
the running time. With the introduction of and ?/8, 
the new definition of vectors uric and vric (see § 12. 2p at 
the star's boundaries becomes 



URi, 



(yi\ 




/2/2\ 


2/4 


, VRic = 


2/3 


2/5 


2/6 


V2/7/ 




\yj 



(62) 



and the new initial conditions on the Riccati matrices 
at the star's center and surface are derived accordingly. 
Even though the introduction of the new variables in- 
creases the size of the various matrices, it does not affect 
the running time significantly. In particular, as far as 
reimbedding f§ I2.2.3P is concerned, during the search for 
the permutation yielding R' with the minimum Eucle- 
dian norm, 1/7 and j/g are kept fixed. This trick yields 
the same number of trials as in the non-adiabatic stellar 
pulsation problem. 
The tidal eigenfunctions are calculated as described in 

§EMl 

3.4. Testing the Extension of the Riccati Method to 
Investigate Dynamic Tides 

Here we test whether our extension of the Riccati 
method to treat tidally excited stellar pulsation works as 
expected. First, we compare the tidal eigenfunctions cal- 
culated with CAFein with the work presented by PS90 in 



the adiabatic regime. Next, we compare the orbital and 
spin evolution timescales due to dynamic tides (§ 13. 2p 
computed wit h CAFein with the results presented by 
iWillems et all (2003) (WVHS03, hereafter). 

3.4.1. Ttdal Eigenfunctions of a 5Mq MS star 

Here, we investigate the effect of non-adiabatic dy- 
namic tides on a MS star of 5Mq (primary) and metal- 
licity Z = 0.018, by studying the variation of the ra- 
dial component of the tidal displacement at the surface 
[^r(l)] as a function of the tidal forcing frequency. Our 
main goal here is to verify numerically that a dynami- 
cal tide can be approximated as the sum of the static 
tide and another part reflecting the oscillatory proper- 
ties of the star itself, in agre ement with what asymp- 
totic theories have shown (e.g. lZahn|[T975l : lSmeverslll997l : 
Emcyers & Willcms 1991). 

As before, we create the stellar model with MESA 
and increase the number of mesh points as described 
in § 12.3.31 The radius of the model adopted here is 
i?i ~ 2.62 _R0, and the convective core extends out 
to r/Ri ~ 0.185. For comparison, these same param- 
eters for the model used by PS90 are Ri ~ 2.52 Rq and 
r/Ri ~0.18, respectively. We set {l,m,k) = (2,0,1), we 
take the companion to be a point mass of 5.0 M©, and 
we set e = 0.4. We then vary the orbital period so that 
the tidal forcing frequency becomes comparable to the 
gg- and gio-mode frequencies. Here we note that since 
TO = and (Tm.fc = fcr^orb + toJIi, the star's spin is not 
relevant. 

The behavior of the radial component of the displace- 
ment from the equilibrium position at the surface |^r(l)l 
is shown in Fig. [5] as a function of the tidal forcing 
frequency (and orbital frequency ). This Figure can be 
compared with Fig. 1 in PS90 or iSavonije fc Papaloizoul 
(1983). The difference bet ween the Figures presented by 
these two studies is that iSavoniie fc Papaloizoul (jl983f ) 
performed a non-adiabatic calculation, while in the adi- 
abatic calculation of PS90 ^r(l) changes discontinuously 
between ± 00. Fig. [5] shows that, close to a resonance, 
the amplitude of the displacement at the surface varies 
greatly for small changes in cjt, while away from a reso- 
nance changes in the amplitude occur at a much slower 
rate. 

In Fig. [6] we show the radial and orthogonal com- 
ponents of the tidal displacement and the total per- 
turbation of the gravitational potential as a function 
of the radial coordinate, for tidal forcing frequencies 
~ 0.4269 - 0.4272. The nine zeros displayed by ^r/r in- 
dicate that we are close to a resonance with the (/g-mode. 
Our Fig. inican be compared with Figs. 2 and 3 in PS90, 
where the discontinuity in ^r(l) is again a result of the 
adiabatic treatment. As PS90 pointed out, the behavior 
of (^r/f )r at the surface is determined by BC (fT8|) and, 
therefore, by the competition of the orthogonal compo- 
nent of the tidal displacement and the total perturbation 
of the potential. This can be seen by considering that, 
close to a resonance with the gg— mode, c± 0.2 and 
uj^'^ ~5.4, while the other term V entering BC (fT8|) is 
y ~ 1.5 X 10^, so V-^ ~ 6.7 X 10-4. Therefore, for 
/ — 2, it can be shown that BC (fTSll reduces to 



(63) 
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Figure 5. Modulus of the radial component of tlie tidal displace- 
ment at the surface for the MS star of 5Mq described in the text. 
The peaks correspond to a resonances with the modes gg and gio. 
For the case under consideration, lot = fcf^orb + mf^i = f^orb a^nd 
the frequencies are in the units of Table [2] 
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Figure 6. Real component s of th e tidal response for the SMq MS 
star described in the text I3.4.ip . ojt is set close to the gg-mode 
eigenfrequency. 

As shown in Fig. [51 close to resonance the behavior 
of ^r(l) is related to the trend of ^h(l), as the contribu- 
tion from the total perturbation of the potential remains 
small. Moving along the steep slope in Fig. [5] towards 
smaller wt , the rapid decrease in ^r(l) is related to the 
rapid decrease of ^h(l) (as PS90 pointed out referring to 
Fig. 1 in their paper). This can be seen from the decrease 
in ^1.(1) in going from = 0.4269 to 0.4261 (from the 
orange line in Fig. [5] to the top plot in Fig. [3 to be 
compared with the behavior presented by PS90 in going 
from their Fig. 3 to 4). At ojt = 0-4261 the nine zeros in 
^i (r)/r are still visible. Along the nearly horizontal line 
of Fig.[5l as cjt decreases further, the increasingly signifi- 
cant total perturbation of the potential affects the behav- 
ior of ^r(l)- As shown in the middle-top panels in Fig. [7] 
(to compare with Fig. 5 in PS90), at uj^ = 0.4058, 6(r)/r 
close to the star's surface shifts from the horizontal dot- 
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Figure 7. Real components o f the tidal response for the 5Mq 
MS star described in the text (§[3ATJ. The horizontal dotted line 
marks the position of the zero, while the dashed line in the left 
panels indicates the static tide component —'^/g divided by r (e.g. 
PS90). In the top and bottom plots, is close to the gg and gio 
eigenfrequency, respectively (see text). 

ted line marking the position of ^i.(r)/r ~ 0. Decreas- 
ing further the tidal forcing frequency to ujt — 0.3972 
(middle-bottom panel in Fig. [7]), the radial component of 
the displacement makes a turn upward at the surface (to 
compare with Fig. 6 in PS90). This behavior is related 
to the change in sign of ^h(l)- Finally, at ujt — 0.3928 
(bottom panel in Fig. [7]), the contribution from the or- 
thogonal component of the displacement and the total 
perturbation of the potential become comparable, and 
the resulting amplitude of ^r(l) is zero (to compare with 
Fig. 7 in PS90). The radial component of the displace- 
ment from equilibrium is approaching the giQ-mode. 

From the behavior of (r) / r displayed in Fig. [3 it is 
clear that away from resonances a dynamical tide can 
be approximated as the sum of the static tide 
denoted with a dashed line) and another part refiect- 
ing the oscillatory properties of the star itself, as a symp- 
totic theories have shown (e.g. iZahn..l975;,Smever3iil997l : 
ISmevers fc Willemsl[T998h . 

3.4.2. Dynamic Tides Timescales in an Eccentric Binary 
Hosting a 5Mq MS star and a Neutron Star 

Here, we test CAFein's results on the orbital and spin 
evolution timescales due to dynamic tides both in and 
out of resonance by reproducing the timescales presented 
by WVHS03 for an eccentric binary hosting a 5 Mq MS 
star (primary) and a neutron star. Note that, for the 
calculation of the dynamic tide timescales, WVHS03 did 
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not solve the fully non-adiabatic problem, as we do here, 
and the perturbed stellar quantities are found via semi- 
analytical solutions (they do use a full non-adiabatic 
calculation for the eigenfrequencies used in the semi- 
analytical solutions). 

As before, we use MESA to create a stellar model of a 
5 Mq MS star at solar metallicity and increase the num- 
ber of mesh points as described in § 12.3.31 The radius 
of the model adopted here is i?i ~ 2.66 i?0, its dynam- 
ical (Tdyn), thermal (rth) and nuclear (rnuci) timescales 
are 323.9 min, 2.71xl0^yr and 9.20xl0^yr, respectively. 
For comparison, these same parameters for the model 
used by WVHS03 are Ri ~ 2.8 i?©, rdy„ =54.8 min, 
Tth =4.88xl0^yr, and Tnuci =8. 67x10'^ yr, respectively. 
In both cases the H fraction at the center is Xc =0.7. 
We take the companion to be a point mass of 1.4 Mq, the 
eccentricity e = 0.5, and we take the spin of the primary 
to be 50% of the companion's orbital angular velocity at 
periastron. 

Similarly to WVHS03, we consider the dominant terms 
in the expansion of the tide-generating potential and fix 
/ = 2 and m — —2. We then calculate the orbital and 
spin evolution timescales for orbital periods (Porb) rang- 
ing from 2 to 5 days, taking into account several terms in 
the expansion of the tide-generating potential and con- 
sidering fc up to 9 in Eq. [51] (WVHS03 considered k up to 
20 and found resonances up to fc = 13). The results are 
summarized in Fig. |8l which should be compared with 
Fig. 3 in WVHS03 (we used the same range for the x- 
and y-axis). It is clear that, even though the resolu- 
tion used by WVHS03 during the scan of the parameter 
space in Porb is higher than the one adopted here, the 
magnitude and trend of the timescales due to dynamic 
tides both in and out of resonance are in good agree- 
ment (and, as pointed out by WVIIS03, also in agree- 
ment with what previous inv estigations have found, e.g. 
iSavoniie fc Papaloizoul [19831) . For the results presented 
by WVHS03, the resonantly excited eigenmodes in the 
range Porb =2-5 days are (7— modes of radial order n 
from 1 to 12. In Table [6] we list the real part of the 
non-adiabatic eigenfrequencies for the modes gi — gi2 of 
the stellar model adopted in this work. As our stellar 
model differs from the one used by WVHS03, a direct 
comparison between their eigenfrequencies (see Table 1 
of WVHS03) and the one calculated here is not possi- 
ble. However, making WVHS03's eigenfrequencies di- 
mensionless, wr for the gi- and gi2-modes calculated by 
WVHS03 are 2.20672 and 0.36723, respectively. This 
suggests that the resonantly excited modes in our stellar 
model are in the range between gi to 511, in agreement 
between the two studies. 

4. SUMMARY, DISCUSSION, AND CONCLUSIONS 

Here we have presented CAFein, a new computational 
tool for calculating non-adiabatic stellar oscillations in 
isolated stars and tidally excited stellar oscillations in 
close binaries, particularly in the dynamic tides regime, 
where the driving frequencies are comparable to the stel- 
lar eigenfrequencies. Tides are considered as a small 
perturbation applied on a spherically symmetric star 
in hydrostatic equilibrium and the linear approxima- 
tion is adopted. CAFein is based on the so-called Ric- 
cati method, a numerical algorithm that has been ex- 



Table 6 

Real component of the non-adiabatic I = 2 
eigenfrequencies for the 5M p, MS described in 

§1332] 



mode 




mode 


t^R 






91 


2.04630 


97 


5.36279 X 


10 


-1 


92 


1.41731 


98 


4.78673 X 


10- 


-1 


93 


1.07746 


99 


4.33268 X 


10- 


-1 


94, 


8.62371 X 10-1 


910 


3.9556 X 


10- 


1 


95 


7.16724 X 10-1 


911 


3.6354 X 


10- 


1 


96 


6.12660 X 10-1 


912 


3.3648 X 


10- 


1 




r 



2.0 2.5 3.0 3.5 4.0 4.5 5.0 
Porb (min) 

Figure 8. Orbital and spin evolution timescales due to dynamic 
tides for a binary hosting a 5 Mq MS star and a neutron star 
in an eccentric orbit (e = 0.5) as a function of the orbital pe- 
riod {Pa^y^).Top: log|ta[ = log|a/ascc|- Middle: logltn^ | = 

log|f2i/r2i^sec|. Bottom: log|te| = log|e/escc|. The horizontal dot- 
ted line represents the logarithm of the star's nuclear time scale. 
A comparison with Fig. 3 in WVHS03 indicates satisfactory agree- 
ment on the timescale calculation. 

tensively and successfully applied to a variety of stel- 
lar pulsators and which does not suffer from the ma- 
jor drawbacks of commonly used shooting and relax- 
ation schemes. Even though the Riccati method is for- 
mally a shooting method, it relies on transforming the 
linear first-order boundary eigenvalue problem describ- 
ing stellar oscillations into a numerically stable, non- 
linear initial value problem. This initial-value problem 
is then solved using a shooting method, where the eigen- 
frequency is the only shooting parameter to be iterated. 

The inclusion of the tide-generating potential in the 
stellar pulsation equations formally does not change the 
system of equations that have to be integrated, and thus 
the applicability of the Riccati method. However, it ren- 
ders the BCs at the star's surface non- homogeneous. We 
made the Riccati method viable for solving the tidally 
excited stellar pulsation problem by introducing two 
new variables and corresponding differential equations 
to make the BCs homogeneous. 

We tested CAFein as a pure stellar pulsation code 
for two different applications. In the adiabatic regime, 
we first calculated the eigenfrequencies of a polytrope 
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and verified that the results are not significantly af- 
fected if some of the relevant parameters entering the 
Riccati method are varied, thus demonstrating CAFein's 
numerical stability. Next, we compared the computed 
eigenfrequencies with previously published results which 
relied both on the Riccati method (TL04) and on 
other commonly used shooting and relaxation techniques 
(CDM94). The comparison yielded very good agreement 
and the orthogonality of low order eigenfunctions was 
also successfully verified. In the non-adiabatic regime, 
we considered a stellar model in the /3 Cephei/SPB insta- 
bility strip of the HR diagram and successfully recovered 
the unstable modes in the part of the parameter space 
examined. 

We showed that the extension of the Riccati method 
to treat tidally excited stellar pulsations works as ex- 
pected by successfully reproducing the work presented 
by PS90 and showing that a dynamical tide can be ap- 
proximated as the sum of the static tide and another part 
reflecting the oscillatory properties of the star itself. Fur- 
thermore, we successfully reproduced the magnitude and 
trend of the orbital and spin evolution timescales due to 
dynamic tides both in and out of resonance presented 
by WVHS03. Here we stress the fact that WVHS03 did 
not solve numerically the fully non-adiabatic problem but 
found the perturbed stellar quantities via semi-analytical 
solutions. 

In this paper we have explored CAFein's performance 
in the dynamic tides regime, where the tidal forcing fre- 
quencies are close to the star's eigenfrequencies and the 
latter can be resonantly excited by tides. This code could 
be used to explore cases that are closer to the quasi-static 
tides limit, in which the tidal forcing frequencies are 
much smaller compared to the inverse of the WD's dy- 
namical time scale (almost synchronized components or 
lo ng orbital and rotationa l perio ds). However, as noted 
bv iSavoniie fc Papaloizoul (|1983[ ). it is notoriously diffi- 
cult to calculate low-frequency tides by integrating the 
full set of tidal oscillation equations because of the short- 
wavelength components entering the tidal response. For 
this reasons, the system of equatio ns is usually re duced 
using perturbation theory (e.g. iWillems et al.l[20rQ' ). 

Even though the physics included in CAFein makes it 
suitable for investigations of stars with envelopes that 
are mostly radiative we intend to use this novel code to 
investigate a variety of binaries and stars. To this pur- 
pose, we are currently upgrading it to account for the 
effect of rotation on the eigenmo de spectruni in th e so- 
called traditional approximation (|Unno et al.l 119890 and 
the effect of turbulent friction acting on the equilib- 
rium tide (iTerauem et allll998t ISavoniie fc Wittell2002l: 
IWillems et al.ll2ninn . 

We are grateful to A. Barker, Y. Lithwick, Lars 
Bildsten, Bill Paxton, and the MESA community for 
useful discussions during the development of CAFein. 
Simulations were performed on the computing cluster 
Fugu available to the Theoretical Astrophysics group at 
Northwestern and partially funded by NSF grant PHY- 
0619274 to VK. This work was supported by NASA 
Award NNX09AJ56G to V.K. 
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